R Markdown

First up is to load in the necessary packages. You may need to install them. brms uses stan which can be a pain to load. It might ask you if you want to compile. If I’m having trouble, sometimes I just say no and it works. Also, start with a clean R environment when you install and load up brms or anything else with stan in it. There are some libraries that can get in the way. I don’t really know the source of the difficulty. I just know sometimes it installs quickly and easily and sometimes it doesn’t. You may need to google around a bit if you’re having trouble.

#install.packages("brms")
#install.packages("BayesFactor")
#install.packages("simstandard")
#install.packages("tidyverse")
library(brms)
## Loading required package: Rcpp
## Loading 'brms' package (version 2.15.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
## 
## Attaching package: 'brms'
## The following object is masked from 'package:stats':
## 
##     ar
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.4     ✓ dplyr   1.0.2
## ✓ tidyr   1.1.2     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(simstandard)
library(BayesFactor)
## Loading required package: coda
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## ************
## Welcome to BayesFactor 0.9.12-4.2. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
## 
## Type BFManual() to open the manual.
## ************

##Data simulation I’ve created four data sets. Two have 40 observations and two have 400. Priors will affect the results differently depending on the number of observations. Each set (small and large) also have a case in which the data will come out to be sig (p<.05) and a case where they will not (p>.05). Note that in the sig case, the data is standardized (the intercept will therefore be close to 0). I’ve posted the histograms so you can see how the distributions look and have computed the correlations between the x (age) and y (vocabulary knowledge) variables to verify that the relationship in the significant set is greater than in the non-sig. set. Using the same random see (set.seed()) will usually ensure you get the same results as I do.

set.seed(4321)
df.ns.small <- tibble::tibble(vk = rnorm(40,200,50), age=rnorm(40,18,4))
summary(df.ns.small)
##        vk             age       
##  Min.   :120.6   Min.   :12.77  
##  1st Qu.:178.1   1st Qu.:17.48  
##  Median :201.1   Median :18.74  
##  Mean   :203.0   Mean   :18.95  
##  3rd Qu.:233.8   3rd Qu.:21.12  
##  Max.   :280.5   Max.   :26.32
hist(df.ns.small$age)

hist(df.ns.small$vk)

cor(df.ns.small$age,df.ns.small$vk) #.02
## [1] 0.01946101
df.s.small <- simstandard::sim_standardized("age ~~ 0.75 * vk", n = 40)
summary(df.s.small)
##       age                 vk          
##  Min.   :-2.00079   Min.   :-1.52273  
##  1st Qu.:-0.81692   1st Qu.:-0.50523  
##  Median : 0.01170   Median : 0.01878  
##  Mean   :-0.06033   Mean   : 0.08363  
##  3rd Qu.: 0.82217   3rd Qu.: 0.82717  
##  Max.   : 2.23053   Max.   : 1.95663
hist(df.s.small$age)

hist(df.s.small$vk)

cor(df.s.small$age,df.s.small$vk) #.64
## [1] 0.6385461
set.seed(4321)
df.ns.large <- tibble::tibble(vk = rnorm(400,200,50), age=rnorm(400,18,4))
summary(df.ns.large)
##        vk             age        
##  Min.   : 77.5   Min.   : 5.898  
##  1st Qu.:170.5   1st Qu.:15.148  
##  Median :202.4   Median :17.909  
##  Mean   :203.0   Mean   :17.804  
##  3rd Qu.:236.2   3rd Qu.:20.478  
##  Max.   :327.0   Max.   :34.077
hist(df.ns.large$age)

hist(df.ns.large$vk)

cor(df.ns.large$age,df.ns.large$vk)#.02
## [1] 0.02486447
df.s.large <- simstandard::sim_standardized("age ~~ 0.75 * vk", n = 400)
summary(df.s.large)
##       age                 vk          
##  Min.   :-2.81939   Min.   :-2.84270  
##  1st Qu.:-0.56397   1st Qu.:-0.63119  
##  Median : 0.02679   Median :-0.02617  
##  Mean   : 0.01714   Mean   : 0.01719  
##  3rd Qu.: 0.66817   3rd Qu.: 0.60855  
##  Max.   : 2.92719   Max.   : 3.52546
hist(df.s.large$age)

hist(df.s.large$vk)

cor(df.s.large$age,df.s.large$vk)#.71
## [1] 0.7136371

##Setting priors I’m going to do everything in brms because I think it’s the easiest to work with all things considered. stan_lm in the rstanarm library is also good.

I’m setting the priors below to be either diffuse (which here is still somewhat informative), weak, or informative based on a normal distribution. If you’re doing logistic regression, you would actually probably want to set the priors on a beta distribution (not the same beta as the beta coefficient) which can be a bit tricky, so I suggest watching a few videos on the beta distribution if you plan to do logistic regression with brms. The cauchy distribution is also commonly used with continuous data. I’m using the normal here because it’s more familiar, but you’ll want to look up the distribution you use to make sure it reflects your prior beliefs.

For the normal, I set the mean and sd. The mean is always the same in this example; it’s the sd that determines how informative the prior is. Class refers to either the variable (b - think beta coefficient), the intercept, sigma (error), or random interecepts and slopes if you’re using mixed models. Coef is the name of the variable (case sensitive). The priors I’m setting up here suggest a belief that for every increase in age by one year, vocabulary knowledge increased by 5 but that the mean (intercept) for vocabulary knowledge was 200 in the case of raw data or 0 in the case of standardized data.

Getting a bayes factor requires setting reasonable priors. The way you set the prior can heavily bias your bayes factor. The best way to handle this is to do a sensitivity analysis to see how your bayes factor reacts to different levels of priors. Here we will also look at how it reacts to priors that suggest there is no relationship (null).

Although we’re primarily interested in the beta coefficient on the variable, setting a prior for the intercept can also be useful. Here I’m using the student_t because it’s normally used for the intercept in default cases. The first number you get is the df, then the mean, then the sd. Because we’re mostly interested in the deviation from the intercept, I think it’s ok to go with the default prior from brms and have only set the prior on the intercept once as an example.

#Priors for an effect, not standardized
diffuse <- c(set_prior("normal(5,100)", class = "b", coef = "age"),
             set_prior("student_t(3,200,2.5)", class="Intercept"))
weak <- set_prior("normal(5,10)", class = "b", coef = "age")
informative <- set_prior("normal(5,1)", class = "b", coef = "age")
hyperinformative <- set_prior("normal(5,.1)", class = "b", coef = "age")

#Priors for an effect, standardized. I only need to do thie for the diffuse prior because it has a prior on the intercept which will change when standardized
diffuse.s <- c(set_prior("normal(5,100)", class = "b", coef = "age"),
             set_prior("student_t(3,0,2.5)", class="Intercept"))

#Null Priors
diffuseN <- c(set_prior("normal(0,100)", class = "b", coef = "age"),
             set_prior("student_t(3,200,2.5)", class="Intercept"))
diffuse.sN <- c(set_prior("normal(0,100)", class = "b", coef = "age"),
             set_prior("student_t(3,0,2.5)", class="Intercept"))
weakN <- set_prior("normal(0,10)", class = "b", coef = "age")
informativeN <- set_prior("normal(0,1)", class = "b", coef = "age")
hyperinformativeN <- set_prior("normal(0,.1)", class = "b", coef = "age")

##Running brms brms is set up just like the lm with a few other parameters. family=gaussian() says that this is a continuous variable and we are using a normal distribution. If this were logistic, I would use family=bernoulli (not binomial). The math for these distributions is different but for inferential purposes, they produce the same results. Brms prefers bernoulli.

Chains refers to the markov chains. 4 is plenty. Explaining this is beyond the scope of this document. Basically, the model starts in 4 places on the distribution and checks to see if it ends up in the same spot each time. rhat=1 in the summary output if this works.

iter refers to the number of iterations in each chain. 2000 is usually enough, but if you find that in the summary output, rhat does not = 1, you should increase your iterations as you don’t have a good fit. rhat = 1.01 is not good enough; you want 1. 2000 is the default, so I have it set once in ns.s.flat so you can see how to do it, but I don’t set it again.

cores allows you to use your computers processing power to speed up the computations. You can use as many cores as there are chains provided your computer has that many. Each core will run a chain. I have two cores on my computer, so if I set cores to 2, the entire process takes ~half the time.

sample_prior = TRUE allows us to plot our priors later. It’s just telling brms not to trash this data. If you’re having trouble understanding your prior, plotting it is the easiest way to get a sense of what is happening. This only works though if you set priors; it doesn’t work with brms default priors.

save_all_pars = save_pars(all=TRUE) keeps the appropriate values on hand for calculating the bayes factor.

I’m running six models below. the first is the basic linear model, then a brms model with flat priors. Flat priors are set by brms. They are uninformative and improper (don’t sum to 1) priors and consist basically of placing equal probability on all real numbers. brms calls this “flat priors over the reals.” We then move on to the diffuse, weakly informative, informative, and hyperinformative models. The next output consists of the prior summaries followed by the results summaries. Results summaries are similar to lm summaries with a few small differences that we’ll go over in class. The final bit are the plots of the brms models. Note that the flat model has no prior distribution plotted.

#Basic linear model
mod.lm.ns.small <- lm(vk~age,data=df.ns.small)
#brms with flat priors (brms sets a prior here that is flat but over the reals meaning only covers the real numbers)
ns.s.flat <- brm(vk~age,data=df.ns.small,family=gaussian(),chains=4,iter=2000,cores=2,sample_prior = TRUE,save_pars=save_pars(all=TRUE))
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
#brms with minimal priors
ns.s.diffuse <- brm(vk~age,data=df.ns.small,family=gaussian(),chains=4,cores=2,prior=diffuse,sample_prior = TRUE,save_pars=save_pars(all=TRUE))
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
#brms with weak priors
ns.s.weak <- brm(vk~age,data=df.ns.small,family=gaussian(),chains=4,cores=2,prior=weak,sample_prior = TRUE,save_pars=save_pars(all=TRUE))
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
#brms with informative priors
ns.s.inf <- brm(vk~age,data=df.ns.small,family=gaussian(),chains=4,cores=2,prior=informative,sample_prior = TRUE,save_pars=save_pars(all=TRUE))
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
#brms with hyper informative priors
ns.s.hyper <- brm(vk~age,data=df.ns.small,family=gaussian(),chains=4,cores=2,prior=hyperinformative,sample_prior = TRUE,save_pars=save_pars(all=TRUE))
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.0/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
##                ^
##                ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.0/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
##          ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1
## Start sampling
prior_summary(ns.s.flat)
##                      prior     class coef group resp dpar nlpar bound
##                     (flat)         b                                 
##                     (flat)         b  age                            
##  student_t(3, 201.1, 44.9) Intercept                                 
##      student_t(3, 0, 44.9)     sigma                                 
##        source
##       default
##  (vectorized)
##       default
##       default
prior_summary(ns.s.diffuse)
##                  prior     class coef group resp dpar nlpar bound  source
##                 (flat)         b                                  default
##          normal(5,100)         b  age                                user
##   student_t(3,200,2.5) Intercept                                     user
##  student_t(3, 0, 44.9)     sigma                                  default
prior_summary(ns.s.weak)
##                      prior     class coef group resp dpar nlpar bound  source
##                     (flat)         b                                  default
##               normal(5,10)         b  age                                user
##  student_t(3, 201.1, 44.9) Intercept                                  default
##      student_t(3, 0, 44.9)     sigma                                  default
prior_summary(ns.s.inf)
##                      prior     class coef group resp dpar nlpar bound  source
##                     (flat)         b                                  default
##                normal(5,1)         b  age                                user
##  student_t(3, 201.1, 44.9) Intercept                                  default
##      student_t(3, 0, 44.9)     sigma                                  default
prior_summary(ns.s.hyper)
##                      prior     class coef group resp dpar nlpar bound  source
##                     (flat)         b                                  default
##               normal(5,.1)         b  age                                user
##  student_t(3, 201.1, 44.9) Intercept                                  default
##      student_t(3, 0, 44.9)     sigma                                  default
summary(mod.lm.ns.small)
## 
## Call:
## lm(formula = vk ~ age, data = df.ns.small)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -83.225 -24.760  -1.386  31.988  77.630 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 198.0078    41.9681   4.718 3.19e-05 ***
## age           0.2625     2.1876   0.120    0.905    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 41.66 on 38 degrees of freedom
## Multiple R-squared:  0.0003787,  Adjusted R-squared:  -0.02593 
## F-statistic: 0.0144 on 1 and 38 DF,  p-value: 0.9051
summary(ns.s.flat)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vk ~ age 
##    Data: df.ns.small (Number of observations: 40) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept   197.71     43.03   112.18   278.62 1.00     3990     2869
## age           0.28      2.23    -4.05     4.78 1.00     3977     2764
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    42.42      4.84    34.27    53.20 1.00     3365     2795
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(ns.s.diffuse)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vk ~ age 
##    Data: df.ns.small (Number of observations: 40) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept   195.90     41.70   115.36   280.19 1.00     3572     2728
## age           0.24      2.19    -4.16     4.49 1.00     3540     2648
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    42.06      4.75    33.76    52.64 1.00     3884     2413
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(ns.s.weak)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vk ~ age 
##    Data: df.ns.small (Number of observations: 40) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept   193.37     41.23   111.02   274.01 1.00     3468     2873
## age           0.51      2.15    -3.70     4.80 1.00     3463     2891
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    42.47      5.00    34.15    53.52 1.00     3253     2802
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(ns.s.inf)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vk ~ age 
##    Data: df.ns.small (Number of observations: 40) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept   123.04     18.65    86.10   158.42 1.00     3653     2697
## age           4.22      0.92     2.43     6.04 1.00     3649     3215
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    43.77      4.96    35.19    54.51 1.00     3347     2708
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(ns.s.hyper)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vk ~ age 
##    Data: df.ns.small (Number of observations: 40) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept   108.40      7.04    94.64   122.18 1.00     3776     3111
## age           4.99      0.10     4.80     5.19 1.00     3664     2631
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    44.36      5.28    35.55    55.95 1.00     3664     2753
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(hypothesis(ns.s.flat, "age = 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Flat Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)

plot(hypothesis(ns.s.diffuse, "age = 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Diffuse Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)

plot(hypothesis(ns.s.weak, "age = 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Weak Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)

plot(hypothesis(ns.s.inf, "age = 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Informative Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)

plot(hypothesis(ns.s.hyper, "age = 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Hyper Informative Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)

The below graphs just put the diffuse, weak, and informative models together into one graph so you can see how the prior distribution affects the posterior.

posterior1 <- posterior_samples(ns.s.diffuse, pars = "b_age")[,c(1,2)]
posterior2 <- posterior_samples(ns.s.weak, pars = "b_age")[,c(1,2)]
posterior3 <- posterior_samples(ns.s.inf, pars = "b_age")[,c(1,2)]

posterior1.2.3 <- bind_rows("prior 1" = pivot_longer(posterior1,c(prior_b_age,b_age)), 
                            "prior 2" = pivot_longer(posterior2,c(prior_b_age,b_age)), 
                            "prior 3" = pivot_longer(posterior3,c(prior_b_age,b_age)), 
                            .id = "id")
modelLME <- lm(vk ~ 1 + age, data = df.ns.small)

ggplot(data    = posterior1.2.3, 
       mapping = aes(x        = value,
                     fill     =  id, 
                     colour   = name,
                     linetype = name, 
                     alpha    = name
                   )) +
  geom_density(size = 1.2)+
  geom_vline(xintercept = summary(modelLME)$coefficients["age", "Estimate"],  #add the frequentist solution too
             size = .8, linetype = 2, col = "black")+ 
  scale_x_continuous(limits = c(-20, 20))+
  coord_cartesian(ylim = c(0, .5))+
  scale_fill_manual(name   = "Densities", 
                    values = c("Yellow","darkred","blue" ), 
                    labels = c("diffuse ~ N(5,100) prior",
                               "weak ~ N(5,10) prior",
                               "informative ~ N(5, 1) prior") )+
  scale_colour_manual(name   = 'Posterior/Prior', 
                      values = c("black","red"), 
                      labels = c("posterior", "prior"))+
  scale_linetype_manual(name   ='Posterior/Prior', 
                        values = c("solid","dotted"), 
                        labels = c("posterior", "prior"))+
  scale_alpha_discrete(name   = 'Posterior/Prior', 
                       range  = c(.7,.3), 
                       labels = c("posterior", "prior"))+
  annotate(geom    = "text", 
           x = 0.45, y = -.13,
           label  = "LME estimate:  0.804", 
           col    = "black", 
           family = theme_get()$text[["family"]], 
           size   = theme_get()$text[["size"]]/3.5, 
           fontface="italic")+
  labs(title    = expression("Influence of Priors on Posterior"),
       subtitle = "3 different densities of priors and posteriors and the LME estimate")+
  theme_bw()
## Warning: Using alpha for a discrete variable is not advised.
## Warning: Removed 3650 rows containing non-finite values (stat_density).

We will do all of the above again, except this time for the data set we know is significant. I’m now using include = FALSE in the header in order to suppress the long progress output within the markdown file.

summary(mod.lm.s.small)
## 
## Call:
## lm(formula = vk ~ age, data = df.s.small)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.76709 -0.37218 -0.07462  0.45164  1.17837 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.1155     0.1060   1.090    0.283    
## age           0.5284     0.1033   5.115 9.27e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6691 on 38 degrees of freedom
## Multiple R-squared:  0.4077, Adjusted R-squared:  0.3922 
## F-statistic: 26.16 on 1 and 38 DF,  p-value: 9.271e-06
summary(s.s.flat)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vk ~ age 
##    Data: df.s.small (Number of observations: 40) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.11      0.11    -0.10     0.33 1.00     3278     2563
## age           0.53      0.11     0.33     0.74 1.00     3162     2944
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.69      0.08     0.55     0.88 1.00     3118     2724
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(s.s.diffuse)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vk ~ age 
##    Data: df.s.small (Number of observations: 40) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.12      0.11    -0.10     0.33 1.00     3365     2981
## age           0.53      0.10     0.32     0.73 1.00     3523     2835
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.69      0.08     0.55     0.88 1.00     3036     2615
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(s.s.weak)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vk ~ age 
##    Data: df.s.small (Number of observations: 40) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.11      0.11    -0.11     0.33 1.00     3941     3054
## age           0.53      0.11     0.32     0.75 1.00     3726     2708
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.69      0.08     0.55     0.88 1.00     3615     2723
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(s.s.inf)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vk ~ age 
##    Data: df.s.small (Number of observations: 40) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.12      0.11    -0.11     0.34 1.00     2868     2509
## age           0.58      0.11     0.37     0.80 1.00     3744     2717
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.70      0.08     0.56     0.88 1.00     3481     3029
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(s.s.hyper)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vk ~ age 
##    Data: df.s.small (Number of observations: 40) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.38      0.70    -1.00     1.77 1.00     3712     2791
## age           4.91      0.10     4.71     5.11 1.00     3254     2745
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     4.60      0.52     3.69     5.79 1.01     3382     2468
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(hypothesis(s.s.flat, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Flat Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)

plot(hypothesis(s.s.diffuse, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Diffuse Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)

plot(hypothesis(s.s.weak, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Weak Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)

plot(hypothesis(s.s.inf, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Informative Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)

plot(hypothesis(s.s.hyper, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Hyper Informative Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)

posterior1 <- posterior_samples(s.s.diffuse, pars = "b_age")[,c(1,2)]
posterior2 <- posterior_samples(s.s.weak, pars = "b_age")[,c(1,2)]
posterior3 <- posterior_samples(s.s.inf, pars = "b_age")[,c(1,2)]

posterior1.2.3 <- bind_rows("prior 1" = pivot_longer(posterior1,c(prior_b_age,b_age)), 
                            "prior 2" = pivot_longer(posterior2,c(prior_b_age,b_age)), 
                            "prior 3" = pivot_longer(posterior3,c(prior_b_age,b_age)), 
                            .id = "id")
modelLME <- lm(vk ~ 1 + age, data = df.s.small)

ggplot(data    = posterior1.2.3, 
       mapping = aes(x        = value,
                     fill     =  id, 
                     colour   = name,
                     linetype = name, 
                     alpha    = name
                   )) +
  geom_density(size = 1.2)+
  geom_vline(xintercept = summary(modelLME)$coefficients["age", "Estimate"],  #add the frequentist solution too
             size = .8, linetype = 2, col = "black")+ 
  scale_x_continuous(limits = c(-20, 20))+
  coord_cartesian(ylim = c(0, .5))+
  scale_fill_manual(name   = "Densities", 
                    values = c("Yellow","darkred","blue" ), 
                    labels = c("diffuse ~ N(5,100) prior",
                               "weak ~ N(5,10) prior",
                               "informative ~ N(5, 1) prior") )+
  scale_colour_manual(name   = 'Posterior/Prior', 
                      values = c("black","red"), 
                      labels = c("posterior", "prior"))+
  scale_linetype_manual(name   ='Posterior/Prior', 
                        values = c("solid","dotted"), 
                        labels = c("posterior", "prior"))+
  scale_alpha_discrete(name   = 'Posterior/Prior', 
                       range  = c(.7,.3), 
                       labels = c("posterior", "prior"))+
  annotate(geom    = "text", 
           x = 0.45, y = -.13,
           label  = "LME estimate:  0.804", 
           col    = "black", 
           family = theme_get()$text[["family"]], 
           size   = theme_get()$text[["size"]]/3.5, 
           fontface="italic")+
  labs(title    = expression("Influence of Priors on Posterior"),
       subtitle = "3 different densities of priors and posteriors and the LME estimate")+
  theme_bw()
## Warning: Using alpha for a discrete variable is not advised.
## Warning: Removed 3670 rows containing non-finite values (stat_density).

We’re repeating the above again for large not significant data set

summary(mod.lm.ns.large)
## 
## Call:
## lm(formula = vk ~ age, data = df.ns.large)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -125.933  -31.983   -0.136   33.682  122.293 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 197.8222    10.7598  18.385   <2e-16 ***
## age           0.2924     0.5893   0.496     0.62    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 47.66 on 398 degrees of freedom
## Multiple R-squared:  0.0006182,  Adjusted R-squared:  -0.001893 
## F-statistic: 0.2462 on 1 and 398 DF,  p-value: 0.62
summary(l.ns.flat)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vk ~ age 
##    Data: df.ns.large (Number of observations: 400) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept   197.59     10.69   176.94   218.56 1.00     4103     3106
## age           0.31      0.59    -0.84     1.44 1.00     4172     3159
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    47.77      1.72    44.46    51.17 1.00     4172     2834
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(l.ns.diffuse)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vk ~ age 
##    Data: df.ns.large (Number of observations: 400) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept   196.38     10.81   175.44   216.96 1.00     4021     2807
## age           0.30      0.60    -0.85     1.48 1.00     4000     2827
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    47.76      1.66    44.65    51.22 1.00     3886     2895
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(l.ns.weak)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vk ~ age 
##    Data: df.ns.large (Number of observations: 400) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept   197.80     10.79   176.18   218.85 1.00     4463     3018
## age           0.29      0.59    -0.88     1.46 1.00     4462     3154
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    47.76      1.71    44.55    51.26 1.00     4132     2803
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(l.ns.inf)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vk ~ age 
##    Data: df.ns.large (Number of observations: 400) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept   176.02      9.40   157.76   194.37 1.00     4511     2829
## age           1.51      0.51     0.49     2.51 1.00     4510     2765
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    47.96      1.76    44.61    51.57 1.00     4301     2784
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(l.ns.hyper)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vk ~ age 
##    Data: df.ns.large (Number of observations: 400) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept   115.91      3.11   109.67   121.96 1.00     3876     2959
## age           4.89      0.10     4.69     5.07 1.00     3978     2932
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    51.21      1.85    47.67    54.95 1.00     3572     3049
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(hypothesis(l.ns.flat, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Flat Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)

plot(hypothesis(l.ns.diffuse, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Diffuse Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)

plot(hypothesis(l.ns.weak, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Weak Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)

plot(hypothesis(l.ns.inf, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Informative Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)

plot(hypothesis(l.ns.hyper, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Hyper Informative Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)

posterior1 <- posterior_samples(l.ns.diffuse, pars = "b_age")[,c(1,2)]
posterior2 <- posterior_samples(l.ns.weak, pars = "b_age")[,c(1,2)]
posterior3 <- posterior_samples(l.ns.inf, pars = "b_age")[,c(1,2)]

posterior1.2.3 <- bind_rows("prior 1" = pivot_longer(posterior1,c(prior_b_age,b_age)), 
                            "prior 2" = pivot_longer(posterior2,c(prior_b_age,b_age)), 
                            "prior 3" = pivot_longer(posterior3,c(prior_b_age,b_age)), 
                            .id = "id")
modelLME <- lm(vk ~ 1 + age, data = df.ns.large)

ggplot(data    = posterior1.2.3, 
       mapping = aes(x        = value,
                     fill     =  id, 
                     colour   = name,
                     linetype = name, 
                     alpha    = name
                   )) +
  geom_density(size = 1.2)+
  geom_vline(xintercept = summary(modelLME)$coefficients["age", "Estimate"],  #add the frequentist solution too
             size = .8, linetype = 2, col = "black")+ 
  scale_x_continuous(limits = c(-20, 20))+
  coord_cartesian(ylim = c(0, 1))+
  scale_fill_manual(name   = "Densities", 
                    values = c("Yellow","darkred","blue" ), 
                    labels = c("diffuse ~ N(5,100) prior",
                               "weak ~ N(5,10) prior",
                               "informative ~ N(5, 1) prior") )+
  scale_colour_manual(name   = 'Posterior/Prior', 
                      values = c("black","red"), 
                      labels = c("posterior", "prior"))+
  scale_linetype_manual(name   ='Posterior/Prior', 
                        values = c("solid","dotted"), 
                        labels = c("posterior", "prior"))+
  scale_alpha_discrete(name   = 'Posterior/Prior', 
                       range  = c(.7,.3), 
                       labels = c("posterior", "prior"))+
  annotate(geom    = "text", 
           x = 0.45, y = -.13,
           label  = "LME estimate:  0.804", 
           col    = "black", 
           family = theme_get()$text[["family"]], 
           size   = theme_get()$text[["size"]]/3.5, 
           fontface="italic")+
  labs(title    = expression("Influence of Priors on Posterior"),
       subtitle = "3 different densities of priors and posteriors and the LME estimate")+
  theme_bw()
## Warning: Using alpha for a discrete variable is not advised.
## Warning: Removed 3653 rows containing non-finite values (stat_density).

and finally again with the large sig data set

summary(mod.lm.s.large)
## 
## Call:
## lm(formula = vk ~ age, data = df.s.large)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.00463 -0.41362  0.00164  0.45531  1.73849 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.005016   0.032358   0.155    0.877    
## age         0.710424   0.034956  20.324   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.647 on 398 degrees of freedom
## Multiple R-squared:  0.5093, Adjusted R-squared:  0.508 
## F-statistic:   413 on 1 and 398 DF,  p-value: < 2.2e-16
summary(l.s.flat)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vk ~ age 
##    Data: df.s.large (Number of observations: 400) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.01      0.03    -0.06     0.07 1.00     3781     2686
## age           0.71      0.04     0.64     0.78 1.00     3719     2595
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.65      0.02     0.61     0.70 1.00     4063     3084
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(l.s.diffuse)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vk ~ age 
##    Data: df.s.large (Number of observations: 400) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.00      0.03    -0.06     0.07 1.00     4230     3332
## age           0.71      0.03     0.64     0.78 1.00     4620     2875
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.65      0.02     0.61     0.69 1.00     3717     2671
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(l.s.weak)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vk ~ age 
##    Data: df.s.large (Number of observations: 400) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.00      0.03    -0.06     0.07 1.00     3898     2346
## age           0.71      0.04     0.64     0.78 1.00     4315     2864
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.65      0.02     0.61     0.69 1.00     4081     3126
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(l.s.inf)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vk ~ age 
##    Data: df.s.large (Number of observations: 400) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.01      0.03    -0.06     0.07 1.00     4351     3140
## age           0.72      0.04     0.65     0.78 1.00     4213     2665
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.65      0.02     0.61     0.70 1.00     3961     2927
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
summary(l.s.hyper)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: vk ~ age 
##    Data: df.s.large (Number of observations: 400) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept    -0.05      0.15    -0.34     0.23 1.00     2428     2131
## age           3.76      0.13     3.51     4.01 1.00     2054     2299
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     2.90      0.15     2.61     3.21 1.00     2042     2295
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(hypothesis(l.s.flat, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Flat Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)

plot(hypothesis(l.s.diffuse, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Diffuse Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)

plot(hypothesis(l.s.weak, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Weak Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)

plot(hypothesis(l.ns.inf, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Informative Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)

plot(hypothesis(l.s.hyper, "age > 0"),plot=FALSE)[[1]] + theme_bw() + labs(title="Hyper Informative Priors",y="Density",x="Beta Distribution") + theme(plot.title=element_text(size=18),axis.title=element_text(size=14),axis.text=element_text(size=12), legend.text=element_text(size=14),legend.title=element_text(size=14))+xlim(-10,10)

posterior1 <- posterior_samples(l.s.diffuse, pars = "b_age")[,c(1,2)]
posterior2 <- posterior_samples(l.s.weak, pars = "b_age")[,c(1,2)]
posterior3 <- posterior_samples(l.s.inf, pars = "b_age")[,c(1,2)]

posterior1.2.3 <- bind_rows("prior 1" = pivot_longer(posterior1,c(prior_b_age,b_age)), 
                            "prior 2" = pivot_longer(posterior2,c(prior_b_age,b_age)), 
                            "prior 3" = pivot_longer(posterior3,c(prior_b_age,b_age)), 
                            .id = "id")
modelLME <- lm(vk ~ 1 + age, data = df.s.large)

ggplot(data    = posterior1.2.3, 
       mapping = aes(x        = value,
                     fill     =  id, 
                     colour   = name,
                     linetype = name, 
                     alpha    = name
                   )) +
  geom_density(size = 1.2)+
  geom_vline(xintercept = summary(modelLME)$coefficients["age", "Estimate"],  #add the frequentist solution too
             size = .8, linetype = 2, col = "black")+ 
  scale_x_continuous(limits = c(-20, 20))+
  coord_cartesian(ylim = c(0, 1))+
  scale_fill_manual(name   = "Densities", 
                    values = c("Yellow","darkred","blue" ), 
                    labels = c("diffuse ~ N(5,100) prior",
                               "weak ~ N(5,10) prior",
                               "informative ~ N(5, 1) prior") )+
  scale_colour_manual(name   = 'Posterior/Prior', 
                      values = c("black","red"), 
                      labels = c("posterior", "prior"))+
  scale_linetype_manual(name   ='Posterior/Prior', 
                        values = c("solid","dotted"), 
                        labels = c("posterior", "prior"))+
  scale_alpha_discrete(name   = 'Posterior/Prior', 
                       range  = c(.7,.3), 
                       labels = c("posterior", "prior"))+
  annotate(geom    = "text", 
           x = 0.45, y = -.13,
           label  = "LME estimate:  0.804", 
           col    = "black", 
           family = theme_get()$text[["family"]], 
           size   = theme_get()$text[["size"]]/3.5, 
           fontface="italic")+
  labs(title    = expression("Influence of Priors on Posterior"),
       subtitle = "3 different densities of priors and posteriors and the LME estimate")+
  theme_bw()
## Warning: Using alpha for a discrete variable is not advised.
## Warning: Removed 3663 rows containing non-finite values (stat_density).

posterior1 <- posterior_samples(ns.s.weak, pars = "b_age")[,c(1,2)]
posterior2 <- posterior_samples(l.ns.weak, pars = "b_age")[,c(1,2)]
posterior3 <- posterior_samples(ns.s.inf, pars = "b_age")[,c(1,2)]
posterior4 <- posterior_samples(l.ns.inf, pars = "b_age")[,c(1,2)]

posterior1.2.3.4 <- bind_rows("prior 1" = pivot_longer(posterior1,c(prior_b_age,b_age)), 
                            "prior 2" = pivot_longer(posterior2,c(prior_b_age,b_age)), 
                            "prior 3" = pivot_longer(posterior3,c(prior_b_age,b_age)),
                            "prior 4" =pivot_longer(posterior4,c(prior_b_age,b_age)),
                            .id = "id")
modelLME.l <- lm(vk ~ 1 + age, data = df.ns.large)
modelLME.s <- lm(vk ~ 1 + age, data = df.ns.small)

ggplot(data    = posterior1.2.3.4, 
       mapping = aes(x        = value,
                     fill     =  id, 
                     colour   = name,
                     linetype = name, 
                     alpha    = name
                   )) +
  geom_density(size = 1.2)+
  geom_vline(xintercept = summary(modelLME.s)$coefficients["age", "Estimate"],  #add the frequentist solution too
             size = .8, linetype = 2, col = "grey")+ 
   geom_vline(xintercept = summary(modelLME.l)$coefficients["age", "Estimate"],  #add the frequentist solution too
             
             size = .8, linetype = 2, col = "black")+ 
  scale_x_continuous(limits = c(-10, 10))+
  coord_cartesian(ylim = c(0, .75))+
  scale_fill_manual(name   = "Densities", 
                    values = c("blue","red","green","orange" ), 
                    labels = c("small.ns.weak ~ N(5,10) prior",
                               "large.ns.weak ~ N(5,10) prior",
                               "small.ns.inf ~ N(5, 1) prior",
                    "large.ns.inf ~ N(5, 1) prior"))+
  scale_colour_manual(name   = 'Posterior/Prior', 
                      values = c("black","purple"), 
                      labels = c("posterior", "prior"))+
  scale_linetype_manual(name   ='Posterior/Prior', 
                        values = c("solid","dotted"), 
                        labels = c("posterior", "prior"))+
  scale_alpha_discrete(name   = 'Posterior/Prior', 
                       range  = c(.7,.3), 
                       labels = c("posterior", "prior"))+
  labs(title    = "Influence of Weak and Informative Priors",
       subtitle = "using a small and large data set")+
  theme_bw()
## Warning: Using alpha for a discrete variable is not advised.
## Warning: Removed 2995 rows containing non-finite values (stat_density).

Below are just null models for each data set. These can be used to get the Bayes factors

Now we’ll compute the actual bayes factors. I’ll walk through two ways to do this but there are others.

The first way involves the library BayesFactor. For this data, we’ll use the function lmBF but you’ll want to explore to see if that is appropriate for your data. lmBF gives us the evidence for the alternative versus a null model. There are other comparisons that can be set up. You’ll need to look into the documentation for your specific needs. A nice feature of lmBF is that it provides an error interval on the BF so you don’t need to run it multiple times to see how reliable it is.

Another way to get a BF is through brms native function bayes_factor. For this, you need to run your null model and then you compare the two. Whichever one goes into the parentheses first is the numerator and the BF you return will be the support for that model. Because the bayes_factor function does not give an error interval, I usually run a few times (here just twice) to ensure the results are equivalent across runs. If they’re not, I increase the iterations in the original model. I then report the average as the BF and the extent of variation as a range or percentage.

lmBF gives the evidence for the alternative. I set bayes_factor to give the evidence for the null. To get values the other way around, just take 1/BF.

#Small n.s. data
bf.s.ns.flat <- lmBF(vk~age,df.ns.small) #.31 (evidence for H1)
## Warning: data coerced from tibble to data frame
bf.s.ns.diffuse <- bayes_factor(null.small.ns.diffuse,ns.s.diffuse) #45.32, 45.93 (evidence for H0)
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
bf.s.ns.weak <- bayes_factor(null.small.ns,ns.s.weak) #5.18, 5.15 (evidence for H0)
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
bf.s.ns.inf <- bayes_factor(null.small.ns,ns.s.inf) #7.23, 7.26 (evidence for H0)
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
bf.s.ns.hyper <- bayes_factor(null.small.ns,ns.s.hyper) #9.56, 9.52 (evidence for H0)
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
#small sig data
bf.s.s.flat <- lmBF(vk~age,df.s.small) #1719.149 (evidence for H1)
## Warning: data coerced from tibble to data frame
bf.s.s.diffuse <- bayes_factor(null.small.s.diffuse,s.s.diffuse) #.04364, .044 (evidence for H0)
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
bf.s.s.weak <- bayes_factor(null.small.s,s.s.weak) #.0048, .00484 (evidence for H0)
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
bf.s.s.inf <- bayes_factor(null.small.s,s.s.inf) #8.67, 8.65 (evidence for H0)
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
bf.s.s.hyper <- bayes_factor(null.small.s,s.s.hyper) #>10000,>10000 (evidence for H0)
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
#large ns data 
bf.l.ns.flat <- lmBF(vk~age,df.ns.large) #.12 (evidence for H1)
## Warning: data coerced from tibble to data frame
bf.l.ns.diffuse <- bayes_factor(null.large.ns.diffuse,l.ns.diffuse) #150.58, 150.17 (evidence for H0)
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
bf.l.ns.weak <- bayes_factor(null.large.ns,l.ns.weak) #16.74, 16.72 (evidence for H0)
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
bf.l.ns.inf <- bayes_factor(null.large.ns,l.ns.inf) #6392.59, 6360.84 (evidence for H0)
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
bf.l.ns.hyper <- bayes_factor(null.large.ns,l.ns.hyper) #>10,000, >10,000 (evidence for H0)
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
#Large s data
bf.l.s.flat <- lmBF(vk~age,df.s.large) #5.85 +- .01% (evidence for H1)
## Warning: data coerced from tibble to data frame
bf.l.s.diffuse <- bayes_factor(null.large.s.diffuse,l.s.diffuse) #<.0001,<.0001 (evidence for H0)
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
bf.l.s.weak <- bayes_factor(null.large.s,l.s.weak) #<.0001,<.0001 (evidence for H0)
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
bf.l.s.inf <- bayes_factor(null.large.s,l.s.inf) #<.0001,<.0001 (evidence for H0)
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
bf.l.s.hyper <- bayes_factor(null.large.s,l.s.hyper) #>10000,>10000 (evidence for H0)
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5

Now that we’ve looked at the bayes factors for these models if our priors suggest an effect, let’s look at them if our priors suggest there’s not an effect. We’re doing this because the BFs were pretty wonky in some spots above because I set a prior that was ridiculous. When set one that is more reasonable, our results are less wonky. I’ll run all the brms in one block here.

plot(hypothesis(ns.s.diffuseN,"age=0"),plot=FALSE)[[1]]+labs(title="Small, NS, Diffuse Prior")+xlim(-1.5,1.5)

plot(hypothesis(ns.s.weakN,"age=0"),plot=FALSE)[[1]]+labs(title="Small, NS, Weak Prior")+xlim(-1.5,1.5)

plot(hypothesis(ns.s.infN,"age=0"),plot=FALSE)[[1]]+labs(title="Small, NS, Inf Prior")+xlim(-1.5,1.5)

plot(hypothesis(ns.s.hyperN,"age=0"),plot=FALSE)[[1]]+labs(title="Small, NS, Hyper Prior")+xlim(-1.5,1.5)

plot(hypothesis(s.s.diffuseN,"age=0"),plot=FALSE)[[1]]+labs(title="Small, Sig, Diffuse Prior")+xlim(-1.5,1.5)

plot(hypothesis(s.s.weakN,"age=0"),plot=FALSE)[[1]]+labs(title="Small, Sig, Weak Prior")+xlim(-1.5,1.5)

plot(hypothesis(s.s.infN,"age=0"),plot=FALSE)[[1]]+labs(title="Small, Sig, Inf Prior")+xlim(-1.5,1.5)

plot(hypothesis(s.s.hyperN,"age=0"),plot=FALSE)[[1]]+labs(title="Small, Sig, Hyper Prior")+xlim(-1.5,1.5)

plot(hypothesis(l.ns.diffuseN,"age=0"),plot=FALSE)[[1]]+labs(title="Large, NS, Diffuse Prior")+xlim(-1.5,1.5)

plot(hypothesis(l.ns.weakN,"age=0"),plot=FALSE)[[1]]+labs(title="Large, NS, Weak Prior")+xlim(-1.5,1.5)

plot(hypothesis(l.ns.infN,"age=0"),plot=FALSE)[[1]]+labs(title="Large, NS, Inf Prior")+xlim(-1.5,1.5)

plot(hypothesis(l.ns.hyperN,"age=0"),plot=FALSE)[[1]]+labs(title="Large, NS, Hyper Prior")+xlim(-1.5,1.5)

plot(hypothesis(l.s.diffuseN,"age=0"),plot=FALSE)[[1]]+labs(title="Large, Sig, Diffuse Prior")+xlim(-1.5,1.5)

plot(hypothesis(l.s.weakN,"age=0"),plot=FALSE)[[1]]+labs(title="Large, Sig, Weak Prior")+xlim(-1.5,1.5)

plot(hypothesis(l.s.infN,"age=0"),plot=FALSE)[[1]]+labs(title="Large, Sig, Inf Prior")+xlim(-1.5,1.5)

plot(hypothesis(l.s.hyperN,"age=0"),plot=FALSE)[[1]]+labs(title="Large, Sig, Hyper Prior")+xlim(-1.5,1.5)

And finally, compute those BF. the BF on the flat prior won’t change, so I skipped that step.

#Small n.s. data
bf.s.ns.diffuseN <- bayes_factor(null.small.ns.diffuse,ns.s.diffuseN) #45.46, 45.60
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
bf.s.ns.weakN <- bayes_factor(null.small.ns,ns.s.weakN) #5.16, 5.12
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
bf.s.ns.infN <- bayes_factor(null.small.ns,ns.s.infN) #1.09, 1.1 
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
bf.s.ns.hyperN <- bayes_factor(null.small.ns,ns.s.hyperN) #1,1.001
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
#small sig data
bf.s.s.diffuseN <- bayes_factor(null.small.s.diffuse,s.s.diffuseN) #.04, .044
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
bf.s.s.weakN <- bayes_factor(null.small.s,s.s.weakN) #.0048, .0044
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
bf.s.s.infN <- bayes_factor(null.small.s,s.s.infN) #.0005, .0005
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
bf.s.s.hyperN <- bayes_factor(null.small.s,s.s.hyperN) #.00009, .03
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
#large ns data 
bf.l.ns.diffuseN <- bayes_factor(null.large.ns.diffuse,l.ns.diffuseN) #150.48, 150.38
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
bf.l.ns.weakN <- bayes_factor(null.large.ns,l.ns.weakN) #14.98, 15.07
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
bf.l.ns.infN <- bayes_factor(null.large.ns,l.ns.infN) #1.75, 1.79
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
bf.l.ns.hyperN <- bayes_factor(null.large.ns,l.ns.hyperN) #1.01, 1.01
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
#Large s data
bf.l.s.diffuseN <- bayes_factor(null.large.s.diffuse,l.s.diffuseN) #<.0001,<.0001
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
bf.l.s.weakN <- bayes_factor(null.large.s,l.s.weakN) #<.0001,<.0001
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
bf.l.s.infN <- bayes_factor(null.large.s,l.s.infN) #<.0001,<.0001
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
bf.l.s.hyperN <- bayes_factor(null.large.s,l.s.hyperN) #<.0001,<.0001
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5